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To analyze phase transitions in a nonequilibrium system we study its grand canonical partition func- 
tion as a lunction of complex fugacity. Real and positive roots of the partition function mark phase 
transitions. This behavior, first found by Yang and Lee under general conditions for equilibrium 
systems, can also be applied to nonequilibrium phase transitions. We consider a one-dimensional 
diffusion model with periodic boundary conditions. Depending on the diffusion rates, we find real 
and positive roots and can distinguish two regions of analyticity, which can identified with two dif- 
ferent phases. In a region of the parameter space both of these phases coexist. The condensation 
point can be computed with high accuracy. 
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The investigation of nonequilibrium systems is a grow- 
ing field in statistical mechanics and currently attracts 
much attention. In this context simple models such as 
driven diffusive systems play a paradigmatic role simi- 
lar to the Ising model in equilibrium statistical mechan- 
ics. These systems establish a simple framework in which 
many phenomena can extensively be studied. Moreover, 
driven diffusive systems can easily be mapped to other 
nonequilibrium models, e.g. of polymer dynamics, inter- 
face growth and traffic flow 0-^]. 

A hallmark of many nonequilibrium systems is the ab- 
sence of detailed balance and the support of stationary 
states with non- vanishing currents. Hence these systems 
build a larger class than respective equilibrium systems 
and phase transitions may appear under less restrict con- 
ditions. For example, it is known that spontaneous sym- 
metry breaking and a first-order phase transition may 
occur in one-dimensional nonequilibrium systems with 
short-range interactions @] . 

In thermal equilibrium the probability measures can 
in principle be expressed through an appropriate ensem- 
ble. For driven systems an equally powerful concept is 
missing. In Ref. || a grand canonical partition function 
for nonequilibrium systems has been introduced for the 
first time. For the definition of this function one uses 
the matrix-product representation (see below) of the sta- 
tionary state. Such a representation is not known for 
every system. But for several models it is known and, for 
instance, its existence is guaranteed for open reaction- 
diffusion models . In this Letter we utilize the partition 
function and show that more concepts from equilibrium 
physics may be applied to nonequilibrium systems. We 
develop a Yang-Lee theory || giving us a very powerful 
method to analyze phase transitions. 

We start from the grand canonical partition function 
Z{x) and study its behavior as a function of complex fu- 
gacity x. Although only real values of the fugacity are 
of physical interest, the analytical behavior of thermody- 
namic functions can be revealed completely only by al- 
lowing the fugacity to be complex. This was first found 



by Yang and Lee for equilibrium systems For finite 
systems the roots of the grand canonical partition func- 
tion Z(x) are in general complex or negative if real. But 
in the thermodynamic limit roots may approach the pos- 
itive real axis. This marks a phase transition; in equilib- 
rium systems the pressure p = ksT limy ^^((l/V) logZ) 
is non-analytical, the density p = (d / '8 'log x){p/k^T) is 
discontinuous, and one can distinguish different phases. 

We show here that the above behavior of the roots 
of the partition function can also be found for nonequi- 
librium systems. By similar reasoning, one can clearly 
define phases and determine their properties and fur- 
thermore transfer other concepts known from equilibrium 
statistical mechanics to nonequilibrium systems. A new 
method to compute the phase transition point with great 
accuracy is established, as well. 

To present our results let us consider an one- 
dimensional stochastic diffusion model with L sites and 
periodic boundary conditions. Each site k may be oc- 
cupied by one particle of type 1 or 2, or may be vacant 
(denoted by 0). Starting from random initial conditions 
with fixed particle densities, p\ and p2, the particles can 
rearrange themselves. The only allowed processes are 
local interchanges conserving the number of particles: 
(a)k(fi)k+i -> (P)k(a)k+i with a,j3 G {0,1,2}. In an 
infinitesimal time interval dr these diffusion processes oc- 
cur with probability g a ^ dr. Here we consider the rates: 

.91,2 = 9, .92,1 = 1, ,9i,o = .9o,2 = f (f) 

All other g a< are zero. The stationary state of this model 
has been studied recently by Monte-Carlo simulations 
and mean-field approximations &. For a snapshot of 
a Monte-Carlo simulation see FigHTl For the parameters 
chosen one observes one macroscopic droplet of particles 
which is surrounded by vacancies with a finite density of 
particles. This is a first sign of a first order phase tran- 
sition. Increasing q the droplet shrinks and disappears 
for q larger than some 9o(p); decreasing q it grows until, 
for q = 1, all particles are bound in the droplet. For 
q > 1 the position of droplet fluctuates and is fixed for 
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FIG. 1. A film of a Monte-Carlo simulation of the model 
for q = 1.2, pi = p2 = 0.2 and L = 100. Particles of type 1 
are grey and those of type 2 black. 



q < 1, i.e. translational invariance is spontaneously bro- 
ken. Similar one- and higher-dimensional models have 
been studied in |10|,p[. 

The probability, V(f3,t), for configurations (3 = 
(Pi, fa, . . . , 0l) at time t follows a Master Equation [Q. 
The stationary probability distribution V s t (/3) can be ex- 
pressed as the trace over algebra elements Dp p|: 



V st ((3) = ±tt(D 01 Dp 2 --.D 0L ) 



The normalization factor, Z, is defined as 

Z= > /• /' ..I) 3r )=tr(C r ') 



J2 tT(D 01 D 02 



(2) 



(3) 



with C = Do + D\ + 1>2- In the following we recog- 
nize that Z plays the role of a grand canonical partition 
function. For Eq. (|^) to hold, the Dp have to fulfill the 
quadratic algebra 



qD 1 D 2 - D 2 D 1 = x 2 D 1 + x x D 2 



DiDq = XiDq 



D D 2 = x 2 D 



(4a) 



(4b) 



(4c) 



The two free parameters x\ and x 2 reflect the freedom to 
choose the two densities of particles. A representation of 
this algebra is known || 

D = Go , Di = xiQi , D 2 = x 2 Q 2 (5) 

with the matrices Q given by 

{Go)ij — SuSij 



Qi = 
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a 3 t 3 

a 4 



, Gi = 
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where we have introduced the notations 

a k = r(2{fc- l} r - {k-2} r ) 

s k t k = {k} r (3r - 1 + (2r - l) 2 {k - l} r 



{k} r 



r k - 1 
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In the following we restrict ourselves to equal densities 
of particles, p\ — p 2 = p, and have xi — x 2 — x. Then 



z = z(x)= Yl x N{f) MG^ 

0u/3 2 ,...,l3 L 



(6) 



takes the form of a grand canonical partition function 
where N((3) counts the number of particles in a configu- 
ration (3 and therefore x can be interpreted as a fugacity. 
The density of particles 1 is given by 



P(x) 



(7) 



For the last expression in Eq. (0) we introduced, analo- 
gously to equilibrium physics, the "pressure" 



P(x) = -logZ(x) 



(8) 



Note, that in this context P is not the physical pressure 
of the particles. The singularities and the analytical be- 
havior of p(x) and P(x) are determined by the roots of 
Z(x). To perform a numerical calculation of tr(C L ) we 
change the ensemble slightly and require that at least 
one vacancy is present. In this case we have to compute 
tr(0oC i_1 ) = (C i_1 )n which is much simpler and one 
has to handle (L/2) x (L/2) matrices only. Accordingly, 
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Z(x) oc IIOe-Xj) 



(9) 



1 in a; with roots Xj. 



is a polynomial of degree L 

To study the analytical behavior of P and p we first 
determine the roots of Z(x) in the complex plane. Since 
the partition function is a real polynomial with positive 
coefficients the roots come in complex conjugated pairs 
and real roots are negative. The roots of Z for two dif- 
ferent choices q > 1 are shown in Fig. ^. For small q the 
roots are lying on an elliptic curve. For q larger than a 
critical value, q CI1 t ~ 2.15, the curve is a hyperbola. 
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FIG. 2. The roots of Z in the complex z-plane. Shown 
are the L — 1 roots for q = 2.0: L = 50 (squares) and L = 100 
(diamonds); q = 2.5: L = 50 (stars) and L = 100 (crosses). 




FIG. 3. The density of roots as a function of the arc length 
for L = 100 and q = 1.2, 1.5, 2.0, 2.5 (from top to bottom). 
The curve of roots in the x-plane is ellipsoidal in the first 
three cases and hyperbolic in the last one. The curves are 
shifted by half of their arclength I. 



Let us treat the elliptic before the hyperbolic case. At 
the end we give some remarks on the case q < 1. 

The elliptic case (1 < q < q cr it)- — In this case, the 
imaginary part of the root closest to the positive real axis 
(e.g. root No. 1 for q = 2.0 in Fig. |J) vanishes for large 
L with L _1 while the real part of this root stays finite 
and reaches a positive value, x(q), in the limit L — > oo. 

To investigate further, let us denote the curve from 
the first to the last root (see Fig. |2|) by c(r) with r G 
[0, 1] where I is the arclength of c. The roots are not 
equally spaced and their density varies with the position 
on the ellipse. Furthermore, with increasing system size 
the density of roots along the curve grows linearly in 
leading order and we can denote the number of roots in 
a line element dr by g(r)LdT. The function gM has a 
non- vanishing large L limit and is shown in Fig. |3|. Since 
in this limit the curve c(r) is as well defined we find the 
pressure (up to an additive constant) 

P(x) = J dr g(r) log(z - c(r)) (10) 

It is a thermodynamic quantity in the sense that its limit 
L —> oo exists. If q < q cr it the curve c closes to give an 
ellipse with a positive intercept of the real axis: c(0) = 
x > and g(0) > in the limit L — > oo. One has to 
distinguish two regions of analyticity, the inner and the 
outer of the ellipse in the complex plane. In both regions 
all physical quantities are continuous and differentiable. 
For a plot of P{x) and p(x) see Fig. ||. 

In the outer region of the ellipse {x > x) the asymptotic 



behavior of P is P{x) = logx and we have p(x) = 1/2 in 
this phase. 

For x in the inner of the ellipse we find a p with p(x) < 
p < 1/2. Crossing the ellipse at x the pressure P(x) is 
not differentiable and p(x) has a finite discontinuity; the 
transition is of first order. The jump in the density is 
related to the density of roots g(0) at the real axis |s) 

1/2 - p = 7T^(0) (11) 

The value of x can easily be obtained by extrapolating 
the real part of the first root for L — > oo. Using Eq. (Q) 
and extrapolating one gets p. 

If the density p is fixed in the interval p < p < 1/2 one 
finds the two phases coexistent in the stationary state; 
the dense one with p = 1/2 (i.e. without vacancies) and 
the other with p = p (see again Fig. |l]). The dense phase 
acts like a bottle-neck and accordingly the current is (q — 
l)p(l-p) = ( ? -l)/4. 

Increasing q to q cr i t , the density of particles above that 
condensation starts, p, reaches 1/2 and the density of 
roots on the positive real axis g(0) vanishes according 
to Eq. (pd|). The phases in the inner and outer region 
become similar and one finds a second order phase tran- 
sition at q = <7 cr it. 

Next we want to fix the density of particles and deter- 
mine the value qo(p) of q below which two phases coexist. 
Therefore, the values of p are plotted versus q in Fig. ||. 
In the same plot we give the mean-field approximation 
for the phase transition q^ F {p) = (1 + 6p)/(l + 2p) as it 
has been calculated in 1. In the latter paper the case 
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FIG. 4. The "pressure" P = (l/L)logZ and the density 
p as a function of the fugacity for q — 1.2 and L = 400. 

p = 0.2 was studied by Monte-Carlo simulations as well 
and the phase transition was found at q^f- c = 1.62 ±0.05. 
With the above method we find the more precise value 
go (0.2) = 1.617 ±0.001. 

The hyperbolic case (q > g C rit)- — Also shown in Fig. ^ 
are the roots for q — 2.5. The curve c(r) is a hyperbola. 
The distribution of roots g(r) (see Fig. ||) is similar to the 
elliptic case. The real and imaginary part of the first and 
last roots, and hence the arclength, I, diverges with L. 
We find only the phase with < p(x) < 1/2 continuous 
and monotonically increasing with x. 

The case q < 1. — In this case the roots lie nearly 
equally spaced on a circle. But the diameter of the circle 
shrinks exponentially with L and in the large L limit all 
roots are located at x — 0. In the limit L — > oo we find 
p(x) — 1/2 for x > 0. Details on the L-dependence and 
geometric properties of c(r) for all above cases will be 
given elsewhere Jl2]| . In physically interesting cases the 
density of particles is fixed in the interval < p < 1/2 
and again one finds again coexistent phases. 

In summary, we have demonstrated that by allowing 
for complex values of the fugacity in the grand canonical 
partition function one can analyze the phases and phase 
transitions in nonequilibrium models. In the model pre- 
sented, a first order phase transition has been investi- 
gated. In particular, the transition point can be com- 
puted with great accuracy. 

Due to the general nature of the method it can also 
be applied numerically or analytically to other models 
]12| where the density of particles controls a phase tran- 
sition, e.g. jamming transitions in traffic related models. 
Furthermore, the method is not constrained to investi- 
gate roots in the fugacity plane. One can study complex 
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FIG. 5. The phase diagram in the p-g-plane. The crosses 
mark the phase transition obtained by the described method, 
the thin line results from mean-field considerations. 

roots of some other parameter (e.g. the input rate in the 
totally asymmetric exclusion process jl3|]) as well ]l2| ]. 

I want to thank M.R. Evans, C. Godreche, T. Hcinzcl, 
D. Mukamel and V. Rittenberg for discussions. 
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